Stability of p-orbital Bose-Einstein condensates in optical checkerboard and square 

lattices 
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We investigate p-orbital Bose-Einstein condensates in both the square and checkerboard lattice 
by numerically solving the Gross-Pitaevskii equation. The periodic potential for the latter lattice 
is taken exactly from the recent experiment [Nature Phys. 7, 147 (2011)]. It is confirmed that the 
staggered orbital-current state is the lowest-energy state in the p band. Our numerical calculation 
further reveals that for both lattices the staggered p-orbital state suffers Landau instability but the 
situation is remarkably different for dynamical instability. A dynamically stable parameter region 
is found for the checkerboard lattice, but not for the square. 
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I. INTRODUCTION 

Orbital physics is important in solid-state material 
due to its key role in understanding many interesting 
phenomena including metal-insulator transition, uncon- 
ventional superconductivity, and colossal magnetoresis- 
tance [lj. However, to fully understand the role of orbital 
degree of freedom in real solid materials is challenging 
because of their complex nature. A quantum degener- 
ate gas in the optical lattice @, Q , which is disorder free 
and highly tunable, is an ideal platform to explore high 
orbital physics as the orbital degree of freedom in such 
an ultracold gas is separate from spin and charge free- 
dom automatically. Moreover, a system of neutral bosons 
loaded into an optical lattice at low enough temperature 
has no counterpart in real quantum materials. Bosonic 
atoms can condensate into non-ground state, opening the 
possibility to explore physics that previously might have 
seemed academic or impossible, e.g., the time-reversal 
symmetry breaking superfluidity in the nodal p band 0- 



Several experimental methods have been developed to 
populate the p and higher orbital bands in optical lat- 
tices 043 • Pioneering experiments have been carried 
out by accelerating the lattice Q, dynamically deform- 
ing the double-well potentials as a single site manipula- 
tion [§4TlT| , and exciting atoms into the higher vibrational 
state along a controlled lattice direction through stimu- 
lated Raman transitions [Hj]. The recent implementa- 
tion of orbital degrees of freedom in checkerboard fl3l — 
UU and hexagonal [l6[ optical lattices truly opens an era 
of exploring orbital phases of quantum matter that have 
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no prior analogues in solid state materials. For exam- 
ple, the experiments of Olschlarger et al [l3|-[l5[ show 
that bosonic atoms are loaded and kept in the excited 
p-orbital bands for nearly as long as the ultracold gases 
can be, thus effectively possessing infinite long lifetime 
in the scale of tunneling. In the experiment, atoms are 
transferred from the s-orbital band to the p-orbital band 
through the changing of the double-well relative depth, 
the time of flight images have illustrated the macroscopic 
occupation in the p-orbital real state and complex state. 
The same group also reported the observation of the ex- 
otic d- and /-band superfluid phases [l4|, EH ■ Dynami- 
cally, fermionic and bosonic atoms are made across from 
lower to high orbital bands in optical honeycomb [17| and 
checkerboard [HI lattices, respectively. 

Although the experiments are done with continuous 
optical potential of periodic oscillation, most of the past 
theoretical work employed the standard tight-binding 
method and approximated the system to a lattice model, 
where the atoms are strictly constrained to the p-orbital. 
Many interesting results have been worked out, such as 
the staggered state as the ground state in thep-orbital 
band [j| , superfluid transition to Mott phase [J [H, EH > 
supersolid quantum pha se [2fj| | and quantum strip order- 
ing in triangular lattices [2l[ (for other interesting studies 
and a brief perspective, see, for example, Ref. |22|). In 
this work, we use the continuous model where little the- 
oretical work has been done besides a variational com- 
putation of the lowest energy in the p-orbital band [23j . 
The continuous model provides a better and complete 
description of the system as it can capture the decay to 
the s band that happens in real experiments and applies 
beyond tight-binding approximation. We focus on the 
stability of the p-orbital state, which is crucial to under- 
stand the p-band superfluidity. To understand the su- 
perfluidity in the s band, stability of the Bloch states in 
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the s band has been analyzed theoretically for different 
lattices — one-dimensional lattices (24|, two-dimensional 
square lattices [25[ , and two-dimensional double- well lat- 
tices [26[ . It was also investigated experimentally for a 
one-dimensional lattice [27| . However, for the p-orbital 
system, to the best of our knowledge, only one paper [28| 
discussed the dynamical instability and it is limited to 
the one-dimensional case. 

In this paper, with the Gross-Pitaevskii (GP) equation, 
we calculate exactly the p-orbital band ground state for 
both the two-dimensional (2D) square optical lattice and 
the checkerboard lattice used in the experiment [l3| . We 
confirm that the lowest-energy state in the p-band is the 
staggered state found in Ref. @ . The Landau instability 
and dynamical instability of this state are investigated. 
For the lattice model approach, such a study of stabil- 
ity is not accurate because the s band is always removed 
from the Hamiltonian to make the condensate impossi- 
ble to decay. Our calculation shows that, for both pe- 
riodic potentials, the staggered state always has Landau 
instability as the state is a local saddle point that can 
decay into the s band. For the dynamical stability, these 
two 2D lattices are very different: our numerical search 
does not find any parameter region for the square lat- 
tice, where the staggered state is dynamically stable; in 
contrast, there exists a parameter region for the checker- 
board potential where the staggered state is dynamically 
stable. This is consistent with the intuitive understand- 
ing that the checkerboard potential offers better stabil- 
ity [29(. That is, on general ground, the checkerboard 
potential may be viewed as a particular configuration of 
the simple double-well lattice potential, and the energy 
gap between the lowest s and the first excited p orbital 
bands is much smaller than that between the p and the 
higher excited bands. Consequently, the first-order decay 
of atoms in the p band due to scattering for the checker- 
board lattice is suppressed by energy conservation ac- 
cording to Fermi golden rule, in contrast with the decay 
for the square lattice of single wells where band spac- 
ings are approximately equal in the tight-binding (i.e., 
the simple harmonic oscillator) limit. 

The paper is organized in the following way. In Sec. 
HI! the general theoretical framework of our calculation 
is given. In Sees. IIIII and IIV1 the results for both the 
square lattice and the checkerboard lattice are presented, 
respectively. Finally, conclusions are drawn in Sec. |Vl 



interaction parameter. The 2D GP equation is 

ihdMr) = \ - ^V 2 + V(r) + g\^\ 2 ] </>(r), (1) 

where V(v + ai) = V(r + a 2 ) = V(r) with ai and a 2 the 
lattice vectors; m is the atomic mass, and g is the interac- 
tion parameter. ?/>(r) is normalized as yj J n \tp(r)\ 2 dr = 1 
where the subscript f2 indicates an integral over the unit 
cell with an area £1 — | a.i 1 1 a.2 1 - 

We are interested in the lowest-energy state in the p- 
orbital band. This type of state must be stationary and 
satisfy the time-independent GP equation 



- ^V 2 + V(r) + g\^\ 2 </>(r) = Mr), (2) 

where /i is the chemical potential. The extended solu- 
tion to the above nonlinear periodic equation has the 
form ip(r) — e ,k r /(r). For the usual Bloch states, / is 
a period function with the same period as that of the 
optical lattice, /(r + ai) = /(r) and /(r + a 2 ) = /(r). 
Besides Bloch states, there are other solutions such as 
the period-doubled solutions [H, |3(J, where / satisfies 
/(r + 2ai) = /(r) and /(r + 2a 2 ) = /(r). For the s band, 
the usual Bloch states always have lower energy than that 
of period-doubled states. But for the p-orbital band, pre- 
vious studies 5, 23] have shown that the period-doubled 
solution has lower energy than the corresponding Bloch 
state due to the extra n phase in p-orbital tunneling. 

For the two types of lattices considered in this work, 
the following two Bloch states are degenerate and have 
the lowest-energy among all the Bloch states, 
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and 



G 



(3) 



(4) 



where G = mbi + nb 2 , bi and b 2 are reciprocal lattice 
vectors, and ki = bi/2 and k 2 = b 2 /2. Substituting 
these two equations into Eq.Q leads to a series of non- 
linear equations of either u or v. We use the subroutine 
f solve of MATLAB to solve these equations. 

There are other types of solutions, which can be sym- 
bolically expressed as either 



II. GENERAL THEORETICAL FRAMEWORK 
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Px±y — /-^ (Px i Py) 



(5) 



We consider the Bose-Einstein condensate of bosons 
in a 2D optical potential with periodicity characterized 
by two lattice vectors to be defined below. To compare 
with realistic three-dimensional experimental systems, 
our model applies to the experiments where a strong trap 
is applied along the third direction. We thus neglect the 
third dimension, which only contributes to the effective 



Px±iy — ^2^ = ' x ^~ Hy) 



(6) 



These solutions are period-doubled states and therefore 
non-Bloch. Without interaction, these non-Bloch states 
would have the same energy as the Bloch states P x and 
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P y . With interaction, these states may break the degen- 
eracy, splitting into either lower or higher energy. To find 
these non-Bloch states, one can similarly substitute Eqs. 
(O and ^ into Eq. © and find the coefficients u and 
v numerically. The coefficients u and v found here are 
in general different from the ones found by substituting 
Eqs. © and g]) into Eq. ©. This is the essential tech- 
nical difference from the work in Ref. [23j. Due to the 
time-reversal symmetry, it is sufficient to consider only 

Px+y and Px-\-iy 

We are interested in the stability of these p-orbital 
states. We know that the lowest-energy state in the 
s band is always stable because it is the lowest-energy 
state of the system. It is imperative to know whether 
the lowest-energy p-orbital state is stable or not. In fact, 
this was already the concern at the beginning of study- 
ing the p-orbital states in cold-atom systems [5j, [6j as the 
decay to the s band seems almost inevitable. However, 
it is possible that interaction may be able to stabilize a 
certain p-orbital state and make it a metastable state. 
The primary purpose of this work is to answer whether 
such a possibility exists. As will be shown later, the state 
Px+iy always has the lowest-energy among all examined 
p-orbital states. Consequently, we will focus on the sta- 
bility of this state. 

To examine the stabilities of a state, one can follow 
the well-known procedure [24j and obtain the Bogoliubov 
equation in momentum q space 



where a z 
have 



is the Pauli matrix. For the state P., 



(7) 



x+iy j 



M = 
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£(q) 
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with 



£(q) 



2m 



iq x ) + (d v + iq y ) 
V(r)- n xy + 2g\P x+iy \ 2 



(9; 



2ai) = 
u q (r 4 



Since P x +i y is period-doubled, we have u q (r 
u q (r + 2a 2 ) = ttq(r) and w q (r + 2ai) = 
2a 2 ) = f q (r). To numerically diagonalize the matrb 
a z M, we expand u and v in Fourier series as it q (r) = 
£ G w G e iGr/2 and u q (r) = £ G v G e iG - r / 2 . The diago- 
nalization of a z M for the phonon modes yields the Bo- 
goliubov excitation of the state P c 



x+iy ■ 



This state hat 



Landau instability if part of its Bogoliubov excitations if 
negative; it has dynamical instab ility if part of its Bo- 
goliubov excitations is imaginary [25]. 



it is described by V(x, y) — Vq cos(x) + cos(y) . For 

this lattice, it is convenient to use the following time 
independent GP equation 



d 2 y ) + V(x, y) + c\^\ 2 ip(x,y) = nip(x,y), (10) 



The above equation has been made dimensionless by scal- 
ing energy with 8E r and length with l/2fc^. In this sec- 
tion, x and y are dimensionless. E r is the recoil energy 
and ItL — 2ir/\ is the wave vector of the laser beam. 
The interaction constant is c = mng/h 2 with m being 
the atom mass, n the BEC density (the average particle 
number per site), and g — 2^/2ttH 2 a s / (am) , where a s is 
the s-wave scattering length and a is the characteristic 
length of the harmonic trap along the z direction. 

Following the procedure described in the above section, 
we have numerically computed three states P x , P x +y , and 



P 



x+iy • 



Fig. [T] illustrates how the energies and chemi- 
cal potentials change with the interaction constant c for 
these three states. It is clear from the figure that state 
Px+iy always has the lowest-energy and the energy gap 
to the other states increases with c. This confirms the 
earlier results obtained with lattice model [5| and varia- 
tional method 23]. It is also shown in Fig. ft] that state 



x+y 



always has the highest energy among the three. 



The phase and density profiles of these three states 
are shown in Fig. [2] These three states not only differ 
in phase but also in density. Since the wave functions 
for both states P x and P x +y are real, their phase can 
only be either or ir. Specifically, the state P x has a 
stripe phase structure while the state P x + y has a square- 
shaped one. The wave function of the state P x +i y is com- 
plex and breaks time-reversal symmetry. Consequently, 
this state has much richer phase structure, which is ev- 
idently shown by the staggered orbital currents in Fig. 
[2je). This feature of staggered orbital currents is the 
most prominent predication in Ref. [f|. 
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FIG. 1: (Color online) (Square lattice) Chemical potential 
(left) and energy per lattice site (right) for the P x state (blue 
solid line), Px+iy state (red dashed one) and P x +y (green 
dashed dot line) in the square lattice. Vo = 0.8E r . 



III. SQUARE LATTICE 

The square lattice can be formed by simply overlapping 
two counter-propagating laser beams. Mathematically, 



As the state P x +i y has the lowest-energy among the 
p-orbital states, we focus on the stability of this state. It 
is examined through its Bogoliubov excitations by diag- 
onalizing the matrix a z M. Our computation finds that 
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FIG. 2: (Color online) (Square lattice) Phase and density 
profile for state P x (a), (b), state P x + y (c), (d), and state 
Px+iy (e), (f), respectively, for the square lattice. Arrows in 
(e) indicate the vortex rotating directions. For phase in (a) 
and (c), the dark blue is 0, the brown is ir. x and y have units 
of l/2k L . V = 0.8E r and c = 0.01. 



the Bogoliubov excitations always have a negative part, 
indicating that the P x +i y has Landau instability and is 
not a metastable state. The situation is more delicate for 
dynamical stability. Figure [3] shows the phase diagram 
of dynamical instability, where the stars mark out the 
region of the (momentum) q space where the Bogoliubov 
excitations are imaginary. It is clear from the figure that 
the stable region of the q space increases as the interac- 
tion c decreases. It is reasonable to expect that the whole 
region be stable when c is small enough. However, within 
our numerical capability, we are not able to identify the 
values of c and Vb for which the state P x +i y is free of 
dynamical instability. As shown in Figs. G^a) and[3Jb), 
even for very small c, there are some regions where the 
excitations are imaginary. This means that state P x +i y 
is dynamically stable only for extremely small values of 
c. We have attempted to calculate the critical c for the 
typical experimental Vo [2| and found that they are of the 
order of 10~ 4 . However, despite of intensive efforts, our 
numerical method is in capable of pinning down the exact 
value of these critical c as indicated by the irregular black 
line in the inset of Fig. where the parameter region 
of dynamical instability for the square lattice is marked 
out. These results imply that it is almost impossible to 
use the square lattice to study p-orbital BEC states ex- 
perimentally as dynamical instability can destroy a BEC 



in tens of milliseconds 27 



IV. CHECKERBOARD POTENTIAL 

The optical lattice used in the experiment is a 
checkerboard potential described by 



V{x,y) 



ee z e 
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— ik^x 



(e z cos(a) + e y 
+ e ie e z (e lkLV 



sin(a)) 



ik L x 



-ik L y\ 



ill) 



where e y and e z are unit vectors in each direction. Here 
x and y are the space coordinates, Ul is the laser wave 
vector, a is the polarization angle to the z direction, e 
is the reflection loss, ij describes the small power dif- 
ference between two interferometers, 9 is the phase dif- 
ference between the beams propagating in the x and y 
directions, and Vq is determined by the laser power. The 
angle a can be used to adjust the degree of anisotropy: 
when a = 7r/5, the energy minimum points of the two p- 
orbital Bloch bands are degenerate. The phase difference 
9 controls the relative depth of a double-well. In the ex- 
periment, bosonic atoms are loaded to the p-orbital band 
by adjusting 6, r) « 0.95, e 0.81, and V = 6.2E r with 
the recoil energy E r = fi 2 k\/2m. 

To have a dimensionless time- independent GP equa- 
tion as in Eq. (fTO)) . we scale energy with AE r and length 
with l/\/2fcj,. In the dimensionless expression, the lat- 
tice vectors of the potential are ai = ^/2n(e x + e y ), 
a.2 = \/2ir(—e x + e y ) and reciprocal vectors are bi = 
(e x + By) j\[2 and \>2 = (—e x + e v )/V%. The quasi- 
momentum and coordinate vector are, k = k x hi + k y \)2 
and r = (irai + ya2)/27r, respectively. 



(a) c=0.001 (b) c=0.005 (c) c=0.01 (d) c=0.05 



0.2 



cT 0.12 



0.04 





/ + * 




+ + 
+ + 




* 




0.12 



0.04 



0.04 0.12 0.2 
<L 



0.04 0.12 0.2 



FIG. 3: (Square lattice) Dynamical stability diagram of state 
Px+iy for the square lattice at Vb = 0.SE r . The region of 
dynamical instability is marked by stars. q x and q y are in 
unites of 2fc_L with a step length 0.02. Since the diagram is 
symmetric with respect to q x and q y , we only draw the upper 
triangle for (a) and (c) and the lower one for (b) and (d). 
The upper limit for q x and q y is 0.25 is due to that Px+iy is 
a period-doubled state. 
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FIG. 4: (Color online) (Checkerboard lattice) Chemical po- 
tential (left) and energy per lattice well (right) for the P x 
state (blue solid line), P x +i y state (red dashed one), and 
P x +y (green dashed dot one) in the checkerboard potential. 
V = 6.2E r . 



We have done the same set of numerical computation 
for the three states P x , P x + y , and P x +i y in the checker- 
board lattice as in the square lattice. Figure [4] shows 
their energies and chemical potentials as a function of 
the interaction constant c. Similar to the square lattice, 



the state P 



x+i y 




FIG. 5: (Color online) (Checkerboard lattice) Phase and den- 
sity profile for state P x (a), (b), state P x + y (c), (d) and state 
Px+iy (e), (f), respectively, in the checkerboard potential. Ar- 
rows in (e) indicate vortex rotating directions. In (a) and (c), 
the dark blue is for phase of and the brown for phase of 7r. 
x, and y have units of \j\plkh- Vo = Q.2E r , c — 0.2. 
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in the checkerboard lattice is found to 
have the lowest-energy too. There is however an evident 
difference, i.e., for the checkerboard, the states P x +i y and 
P x are very close in energy while state P x + y has much 
higher energy. This feature suggests that the state ob- 
served in the experiment [l3| is probably not P x + y . 

The phase and density profiles of these three states are 
illustrated in Fig. [5j The phase profiles show a structure 
similar to that for the square lattice, such that P x has 
the wavelike stripe phase structure, P x + y has the square 
phase structure and P x +i y has the staggered orbital cur- 
rents structure This result supports the conclusion 
that one may use the square lattice as a simplified theo- 
retical model to understand much of the unconventional 
properties of the p-orbital condensates as observed in the 
more complex, experimentally realized checkerboard lat- 
tice. In terms of dynamical instability, the two lattice 
configurations are however qualitatively different, to be 
elaborated below. The density profiles do not show clear 
difference between state P x and P x +i y and the reason is 
that the probability density in the deeper well where the 
vortex appears is very small. 

For stability, we focus on that of state P x +iy just as 
in the square lattice. We investigate it also through Bo- 
goliubov excitations by diagonalizing the matrix a z M. 
Similar to the square lattice, our calculation shows that 
the Bogoliubov excitations in the checkerboard always 
have negative part, indicating that state P x +i y in the ex- 
periment also has Landau instability. For dynamical sta- 
bility, the unstable region in the q space increases with 
c as shown in Fig. [5] However, there is a crucial differ- 
ence from the case of square lattice: there exists a critical 
value of c, below which there is no dynamical instability 
as indicated by the blank panel in Fig. [HJa). We are 
able to mark out a region in the space spanned by the 

system parameters c and Vo/E r , where the P x +i y state is free of dynamical instability (shown in the phase diagram 
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FIG. 6: (Checkerboard lattice) Dynamical stability diagram 
of the state P x +i y of the checkerboard potential. The region 
of dynamical instability is marked by stars. The blank panel 
in (a) indicates that the system is free of dynamical instability 
for c = 0.1, to be contracted with that for the square lattice in 
Fig. [3] q x and q y are in unites of V2kL with a step length 0.02. 
Vo = 6.2E r . Since the diagram is symmetric with respect to 
q x and q y , we only draw the upper triangle for (a) and (c) 
and the lower one for (b) and (d) . The upper limit for q x and 
q y is 0.25 is due to that P x +i y is a period-doubled state. 
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Fig. [7]). Due to the uncertainty of the BEC density, we 
have marked the experimental parameter range [l3j with 
a solid line. This shows that it is likely that the P x +i y 
for the experimental setup is dynamically stable. Since 
the time scale for Landau instability is of the order of 
500 m s 12711 which is much longer than that of the exper- 
iment [13], it is reasonable that Landau instability does 
not have much effect. 

In order to make sure that our calculation is correct, we 
simulate the real system of Eq.([T|) with the split-operator 
method. We evolve numerically a BEC in the P x +i y state 
with a small perturbation dtp (10%). When c = 0.2, the 
simulation shows that the state is stable. When c = 
3.0 and c = 7.9, the simulation shows that the state is 
destroyed after t = 17.5ms and t — 2.5ms, respectively. 
All results in the three simulations are consistent with 
our Bogoliubov excitation calculation. 

To map out the phase diagram of dynamical instability 
in Fig. [7] experimentally, one may need to use the Fesh- 
bach resonance to tune the interaction strength. When 
the Feshbach resonance is not available, one can still ob- 
serve the effect of dynamical instability by turning up the 
laser power to drive the system into dynamically unsta- 
ble regime. The effects of dynamical instability should 
be similar to what was observed in Ref. 12711. 
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FIG. 7: (Color online) (Checkerboard lattice) Stable and un- 
stable region in the system parameter space for state P x +iy in 
the checkerboard potential. The black line indicates the pa- 
rameter range used in the experiment [l3| where Vo = 6.2E r . 
The inset shows the stability regions for the square lattice. 
The irregularity of the solid line in the inset is caused by 
the inability of our numerical method to compute precisely 
the critical value of c below which the system is dynamically 
stable. 



CONCLUSION 



In conclusion, we have examined a cold gas of interact- 
ing bosonic atoms loaded in two optical lattice geome- 
tries, namely, the square and the checkerboard lattices, 
for which unconventional p-orbital Bose-Einstein conden- 
sates have been under active investigation in recent years, 
both theoretically and experimentally. The usual theo- 
retical approach used in the past is to assume the stan- 
dard tight-binding approximation and conveniently re- 
duce the system to a Hubbard-like lattice model of one 
single orbital band of interest, i.e., the p band. The 
present approach is however different. Here, the model 
system is solved numerically with the GP equation of 
microscopic two-body interaction by treating the opti- 
cal lattice exactly as a continuous, periodic potential, in 
which both the ground state s-band and all the higher 
orbital bands are included. The approach thus is capa- 
ble of providing a complete analysis for the Landau and 
dynamical instabilities of a p-orbital BEC. Such a com- 
plete analysis for instability was not considered before, 
to the best of our knowledge. We find that the staggered 
state P x +iy indeed has the lowest-energy in the p band. 
By computing the Bogoliubov excitation, we further find 
that for both lattices Landau instability is present, which 
shows that the staggered state is not really a state at lo- 
cal energy minimum. For dynamical stability, we find 
that there exists a parameter region where the staggered 
state is free of dynamical instability for the checkerboard 
lattice whereas no such a parameter region is found for 
the square lattice. This suggests that the staggered state 
be of long life time in the former, but not the latter. 
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